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Abstract: 



A variational approach, based on a discrete representation of the chain, is used to calculate free 
energy and conformational properties in polyelectrolytes. The true bond and Coulomb potentials are 
' approximated by a trial isotropic harmonic energy containing force constants between all monomer- 

^ . pairs as variational parameters. By a judicious choice of representation and the use of incremental 

matrix inversion, an efficient and fast-convergent iterative algorithm is constructed, that optimizes 
the free energy. The computational demand scales as rather than as expected in a more naive 
approach. The method has the additional advantage that in contrast to Monte Carlo calculations 
the entropy is easily computed. An analysis of the high and low temperature limits is given. Also, 
the variational formulation is shown to respect the appropriate virial identities. The accuracy 
of the approximations introduced are tested against Monte Carlo simulations for problem sizes 
ranging from iV = 20 to 1024. Very good accuracy is obtained for chains with unscreened Coulomb 
interactions. The addition of salt is described through a screened Coulomb interaction, for which the 
accuracy in a certain parameter range turns out to be inferior to the unscreened case. The reason 
is that the harmonic variational Ansatz becomes less efficient with shorter range interactions. 

As a by-product a very efficient Monte Carlo algorithm was developed for comparisons, providing 
high statistics data for very large sizes - 2048 monomers. The Monte Carlo results are also used 
to examine scaling properties, based on low-T approximations to end-end and monomer- monomer 
separations. It is argued that the former increases faster than linearly with the number of bonds. 
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1 Introduction 



Polymers and polymer solutions play a profound role in our daily lifc, biologically and technolog- 
ically. This is of course one reason for the intense theoretical studies of polymers, but they have 
also because of their very nature been a challenge to theoreticians. Much theoretical work has been 
done in order to obtain a general understanding of neutral polymers, either in melts or in solution 
10, H. Polyelectrolytes, on the other hand, have been less thoroughly investigated, despite their 
importance in many technical applications - for example in glue production or in food industry, for 
promoting flocculation or in pulp drying ^, ^. The term polyelectrolyte is sometimes used as a 
collective name for any highly charged aggregate. However, here we will restrict the meaning to a 
flexible molecule with several or many charged or chargeable sites; poly-L-glutamic acid, polyamines 
and polysaccharides are some typical representatives. 

The conformation of a flexible polyelectrolyte is a result of the competition between the covalent 
bonding forces, electrostatic interactions as well as more specific short ranged interactions. For 
example, poly-L-glutamic acid and several polysaccharides undergo helix to coil transition as a 
function of pH H, Q . This transition is obviously governed by electrostatic forces - similar structural 
transitions are seen for DNA ||] . Undoubtedly, both the polymer nature and the interaction between 
charged amino acids play an important role for the folding of a protein, as well as other solvent 
averaged forces. The present study should be seen as an attempt to approach the folding problem 
using what turns out to be a very powerful statistical mechanical variational technique. In the first 
step we will limit the study to linear polyelectrolytes in salt solution. 

Variational methods are standard techniques in quantum mechanics, but less so in statistical me- 
chanics, although variational principles were formulated many years ago One type of variational 
formulations starts off from an approximate free energy, which is optimized with respect to the par- 
ticle density [ p^ . For polymers a more fundamental approach is possible, by introducing variational 
parameters directly into the appropriate Hamiltonian. In the past this route has been followed in a 
number of polymer studies |ll|, ^ and it has recently been revitalized by several groups |l4| |l5| . 
To the best of our knowledge, all these calculations have been concerned with continous chains and 
only one or at most a few variational parameters have been optimized. 

The present approach, of which some results were already published |l6), is inspired by refs. ^ 
[l9| . It uses a discrete representation of the polymer, which not only allows us to investigate linear 
or cyclic polymers, but also polymers of arbitrary topology. Thus, e.g., a hyperbranched dendrite 
structure can easily be handled within this formalism. It relies on a variational Ansatz in the form 
of a generic Gaussian distribution, with adjustable force constants between every pair of monomers 
(not to be confused with a Brownian model with forces only between nearest neighbours). Thus, the 
number of variational parameters is proportional to A'^^, where N is the number of interacting units. 
In general, the variational approach is expected to be most accurate at high dimensions ||l^, |l8|. 
Apparantly, this has discouraged the community from pushing the approach for three-dimensional 
polymers into a numerical confrontation. Also, using the method in a naive way would give a 
computational effort scaling like A"^, which would make the method less tractable for large sizes. 
We have also found empirically that such a naive implementation is plagued with bad convergence 
properties. In previous work [ p^ an algorithm was developed that lowers the computational costs 
to A^^ with controlled and nice convergence properties. 

The high and low T limits of the variational approach, which are accessible with analytical means, 
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are derived in this paper. The low T expansions for various quantities yield results that very well 
approximate what emerges from the corresponding expansions in the exact theory. The fact that the 
temperature of interest (room temperature) is fairly low on the temperature scale partly explains 
the success of initial numerical explorations of the variational approach ]l6| and motivates further 
studies. 

Besides being more realistic, the discrete chain also offers the possibility of direct comparison with 
Monte Carlo (MC) simulations, thereby giving an important indication of the accuracy. The final 
output of the present variational calculation is of course the minimized free energy, but one also 
obtains a matrix with all possible monomer-monomer correlations within the chain. This matrix is 
the starting point for the calculation of end-end and monomer-monomer separations and different 
kinds of angular correlations. 

MC simulations of flexible polyelcctrolytes have only recently appeared in the literature and then 
limited to rather short chains [^0[ |2^, |2^, |2^, |2^, , the longest chains being of the order 
of a few hundred monomers. So far simulations have dealt with both explicit representation of salt 
particles as well as implicit in the form of a screened Coulomb potential. The inclusion of short 
ranged interactions has also been studied as have the conformational changes associated with the 
titration of a flexible polyelectrolyte ||2^ . Most simulations have been carried out in the canonical 
ensemble, although the latter problem required the use of grand canonical MC simulations pTf . 
The accuracy of the screened Coulomb potential has been investigated by several people and found 
to be an excellent approximation for not too high polyelectrolyte charge densities and in the absence 
of any multivalent ions ^ . In order to obtain high statistics MC results a pivot algorithm and 
a recently developed hybrid scheme |]30| were used. As a by-product these simulation results were 
also used to examine certain scaling properties, theoretically derived from an approximate low T 
expression. In contrast to the case of rigid bonds we find that the end-end separations scale faster 
than linearly with N. 

When confronting the variational approach with MC data the following results emerge in this work: 



• The variational approach has the unique property that it directly yields the free energy F . 
This is in contrast to MC simulations, where F is only indirectly accessible through elaborate 
integrations. For N=20, where comparisons are inexpensive, the variational free energy nicely 
agrees with MC results. 

• For unscreened Coulomb chains the success reported in ref. [|l6) survives to even larger systems 
(N=1024) for configurational properties like end-end correlations. Also for angular correlations 
and scaling behaviour the variational approach reproduces MC data very well. 

• When including salt through Debye screened Coulomb potentials the performance of the 
method deteriorates somewhat on the quantitative level. Even if the qualitative configura- 
tional picture agrees with that from MC data, the actual numbers for e.g. end-end correlations 
could differ up to 50 %. We attribute this discrepancy partly to the inability of harmonic forces 
to reproduce short range interactions. 



One should also mention that other approximation schemes than the variational ones have been 
attempted in order to study polyelectrolyte conformations. In the mean field approximation it is 
only with the assumption of spherical symmetry that the mean field equations become tractable. 



2 



but the results are not particularly encouraging In a cylindrical geometry, however, the mean 
field approximation behaves much more satisfactorily, but at the expense of large numerical efforts. 

This paper is organized as follows: In sect. 2 the basic Coulomb models (screened and unscreened) 
are presented. The variational approach is presented in sect. 3. The high and low T limits of 
the variational scheme are computed with analytical methods in sect. 4. A description of the MC 
method used in order to establish the quality of the variational method is found in Sect. 5. The 
results from confronting the variational approach with MC data are presented in sect. 6. Finally 
in sect. 7 a brief summary and outlook is given. Most of the detailed derivations are found in 
appendices - generics about the variational method (A), variational energies for a polyelectrolyte 
(B), virial identities (C), high and low T expansions (D) and zero temperature scaling properties 
(E). The disposition of the material into bulk text and appendices is such that the approach and 
results can be fully understood without reading the appendices. 



2 The Model 



2.1 The Unscreened Coulomb Chain 



In this model we consider a polyelectrolyte at infinite dilution and without any added salt. The 
polyelectrolyte counterions are thus neglected and the only electrostatic interactions are between 
the charged monomers. More explicitly, the polymer chain consists of N point charges connected 
by harmonic oscillator ("Gaussian") bonds. The potential energy for a chain then takes the form 

Here, the position of the ith charge, and 

Xij = Xi — Xj (2) 

while q is the monomer charge and e^eo is the dielectric permittivity of the medium. We use the 
tilde notation E, Xj, etc. for physical quantities in conventional units, and reserve E, x;, etc. for 
dimensionless ones, which will be used in the theoretical formalism below. 

The force constant can be reexpressed in terms of the N — 2 equilibrium distance tq, given by 

0^ 1 

k = 3 (3) 

In terms of the dimensionless coordinates x^ , defined by 

Xi = roXi (4) 

the energy takes the form 

E^krt I 1 V|x„+ir + V — I (5) 
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T = ^ (6) 
kri 



We will consider the system at a finite temperature T, which can be similarly rescaled, 

ksT 

with fcs the Boltzmann constant. One then obtains for the Boltzmann exponent the simple expres- 
sion 

E 



In other words, a polyelectrolyte at infinite dilution represents a two-parameter model where T and 
the number of monomers N are the only two non-trivial parameters of the system. 

Unless otherwise stated the following parameter values will be used throughout the paper; T=298K, 
6^=78.3 and ro=6rA. Obviously, E depends only on the relative positions; the global center-of-mass 
position variable will have to be excluded from integrations over the coordinate space. 

The Gaussian and Coulomb energies are subject to a virial identity (see Appendix C), which in 
dimensionless units reads 

2{Eg) ~ (Ec) ^ 3iN - 1)T (8) 

Eq. (||) is a useful relation for checking the correctness and convergence behaviour in MC simulations 
and we find that it is in general obeyed to 0.3% or better. 



2.2 The Screened Coulomb Chain 



To treat a single polyelectrolyte in a solution at finite salt concentration becomes very costly in 
a MC simulation, since for reasonable salt concentrations the number of salt ions easily becomes 
prohibitively large, much larger than the number of polyelectrolyte monomers. The usual way to 
avoid this problem is to preaverage the degrees of freedom of the simple salt ions for some fixed 
configuration of the polyelectrolyte , thus defining salt averaged effective potentials - this is the 



basis of the electrostatic contribution in the classical DLVO potential |32 . In this way we may 
derive a screened Coulomb potential from a linearisation of the Poisson-Boltzmann equation for the 
salt. Eq. (|l|) is then replaced by 



where k is the Debye screening length for a 1:1 salt defined as 



i^ = i\ — (10) 



In eq. (10) Cg is the salt concentration in molars (M) and Na is the Avogadro's number. The 
Boltzmann factor will for the screened Coulomb potential contain the inverse dimensionless Debye 
screening length, k — rgk, as an additional parameter, and with the parameter values given above 
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we have for Cs = 0.01 M, 0.1 M and 1.0 M the ^-values 0.1992, 0.6300 and 1.992 respectively. Then 
eq. is modified into 

E _ 1 
T ~ T 

The virial identity (see Appendix C) now takes the form 

2{Eg) - (Ec) - e^'^''^'^') = HN - l)T (12) 



2.3 Relative coordinates 

In the remainder of this paper, relative coordinates will be used; instead of the absolute monomer 
positions Xj, the bond vectors r.^, 

ri=x.,+i-Xi, i = l,...,7V-l (13) 

will be used as the fundamental variables. In this way complications due to the translational zero- 
mode are avoided; in addition the convergence of the algorithm is considerably speeded up, especially 
at high temperatures. 

The energy of the screened Coulomb chain will then take the following form: 

N — l 

E{v)^Eg + Ec^-Y.v^^+Y. (14) 

where o runs over contiguous non-nil sub-chains, with 

corresponding to the distance vector between the endpoints of the subchain. The unscreened chain 
results for k = 0. 



3 The Variational Approach 

3.1 The Gaussian Ansatz 

In refs. ||l^, ^ |l^ the variational method of refs. [|[ |l^ (see Appendix A for a generic description) 
was revisited in the context of discrete chains of polyelectrolytes. The approach is based on an 
effective energy Ansatz Ey^ given by 

^^/^ = ^ E G - ■ (^^- - (16) 



te- -g^^ 
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where defines an average bond vector, around which Gaussian fluctuations are given by the 
symmetric positive-definite correlation matrix dj , the matrix inverse of which appears in the energy. 

Using this effective energy, the exact free energy F = — Tlog Z of the polymer is approximated from 
above |^ by the variational one 

F = Fv + {E- Ev)v > F (17) 

where Fy = — TlogZy, and ()y refers to averages with respect to the trial Boltzmann distribution 
eM-Ev/T). 

The parameters Gij and are to be determined such that the variational free energy F is minimized; 
we note that the number of variational parameters increases with N like N'^ . The resulting effective 
Boltzmann distribution is then used to approximate expectation values (/) by effective ones (/)y. 
Thus, we have e.g. 

{r^)v = a, (18) 
(r^ • rj)v = Sii ■ SLj + 3Gij (19) 

For the polyelectrolyte systems treated in this study we will at high temperatures find a unique 
variational solution, characterized by a^ = 0; this defines a purely fluctuating solution. At low tem- 
peratures we find in addition a rigid solution with aligned a^ 7^ 0. The latter is due to spontaneous 
symmetry-breaking; it ceases to exist at high temperatures, but will at low enough temperatures 
have the lower free energy. As discussed below, the rigid solution can typically be disregarded at 
normal temperatures. 

For potentials more singular than 1/r^, {E)v will be divergent, and the approach breaks down. 
However, such potentials are not physical and we do not consider this limitation of the approach a 
serious one. 

A non-trivial result of the scaling properties of the effective energy is that the virial identity, eqs. 
(^,^2|), will be respected by the above variational approach (see Appendix C). 



3.2 Using Local Fluctuation Amplitudes 

The minimization of F with respect to Gy and a^ gives rise to a set of matrix equations to be solved 
iteratively. These are considerably simplified, and the symmetry and positivity constraints on Gij 
are automatic, if Gij is expressed as the product of a matrix and its transpose: 

N-l 

Gij — ^ ^ ^ifi^jfi — ^i ' (20) 

The interpretation of the local parameter Zi is simple - it is a fluctuation amplitude for the zth bond 
vector r.i. We can write 

Ti = SLi +^Zi^Jf, (21) 
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where each component of G TZ^ is an independent Gaussian noise variable of unit variance. 
Similarly, we have for a subchain 

r<T = ^r.; = + ^Zct^J^, (22) 

where a.^ = J2iecr ^« ^^"^ ^o- = J2iea Thus, the noise amplitudes are additive. 
The matrix inverse of G can similarly be decomposed: 

G-] = • w, (23) 
where is the (transposed) matrix inverse of Zi^: 

Zi ■ Wj = 6ij (24) 
Note that z^, and are vectors, not in R^, but in R^~^. 

The equations for a local extremum of F(a., z) are obtained by differentiation with respect to and 

OF dP 

^ = 0,^=0 (25) 
3.3 The Unscreened Coulomb Chain 

In terms of and z,, the variational free energy for the pure Coulomb chain becomes, ignoring 
trivial additive constants (see Appendix B): 



F = -3riogdet z + \ Y^{2,zl + a^) + ^ 1 erf 



(26) 



The equations for a minimum will be 

(28) 



dF _ _ \p ^ 



-— - erf I '^'^ 

TT \V2 Za 



where the reciprocal vector is defined by eq. (|24|). 

These equations allow a purely fluctuating solution with a^ = 0. Setting a,; = the variational free 
energy simplifies to 

F = -3Tlogdetz+^^z2 + ,/l^l (29) 

Z . ^ TT Zfj 

% a 

which looks very much like the energy of an (iV — l)-dimensional Coulomb chain with bonds z^, but 
with an extra entropy term (the first) preventing alignment of the ground state. The z derivatives 
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become 



d^i V TT ^ Z% 



(30) 



In order to bring out the structure of the variational solution, we take the scalar product of eq. ( po[ ) 
with Wj to get an expression for the force constants, 

TG^'-S.-l\llT.^'^ (31) 
Thus the variational energy that minimizes the free energy has the following structure: 

The first term is just the original bond term, while in the second term the Coulomb interactions are 
replaced by repulsive harmonic forces, having the right scale to give a good approximation to the 
Coulomb interactions for typical distances Va- (x Za- 



3.4 The Screened Coulomb Chain 



In the case of Debye screening the expression for F is modified to 



3Tlogdetz+i^(3z,2 + a2) (33) 



F „ , 

2 



2 
r 

V / 1 V Z(j j V Zn 



where 

*(x) = exp(xV2) erfc(x/\/2) (34) 
The corresponding derivatives (cf. eq. (|2^)) take the form 

dP 

—— = -3Twi + 3zi (35) 



OF 



Za J V TT KZa 



(36) 



Also here, a purely fluctuating solution is allowed. Setting = the variational free energy reduces 
to 



F = -3Tlogdet + ^ 51 ^» + H I \/f ^ " (''■"-)| 



(37) 
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The z derivatives will be 



g = _3Tw,; + 3z,; - ^ J ^^(1 - n'zl) + K^zl^ (kz,) = (38) 



3.5 Implementation 



Due to the use of relative coordinates and of local noise amplitudes, a simple gradient descent 
method with a large step-size e can be used, that gives fast convergence to a solution of eqs. (ESh 



OF OF 
Az, ^-e,j—, Aa, ^ -eaj—, (39) 

Further speed is gained by updating the reciprocal variables using incremental matrix inversion 
- the increment in Wj due to Az^ is given (exactly) by 

w,-(wi • Az,-) 

1 + w, -Az, ^ ' 

to be applied in parallel for j for fixed i. As a by-product the denominator (1 + w.^ • Az^) gives the 
multiplicative change in the determinant det z, needed to keep track of F. 



In order to maintain a reasonable numerical precision, the erf-related functions needed in the process 
are evaluated using carefully defined Taylor expansions or asymptotic expansions, depending on the 
size of the argument. 



As discussed above, the equations for a minimum are consistent with a purely fluctuating solution 
= 0. Such a solution does indeed exist at all T; furthermore it is the only solution for high 
enough T. It turns out that for realistic choices of T one is in the region where this solution gives 
rise to good results (see figs. |^ below). The additional solution with a^ 7^ appearing at low T is 
a symmetry-broken solution; to be realistic, such a solution should show an anisotropy also in the 
fluctuations, i.e. different amplitudes z^ for the fluctuations parallel and transverse to the direction 
defined by (the aligned) a^. This requires a more general Ansatz than eq. (|l^); theoretically, this 
will produce better low-T solutions, but for the Coulomb chain it leads to equations containing 
functions that arc difficult to evaluate numerically, so we will not use this possibility in this paper. 
The incomplete symmetry-breaking partly explains the tendency for the a 7^ solutions to produce 
inferior solutions. 



For the above reasons, and for reasons of continuity, we will in the numerical explorations use the 
a = solutions, where not otherwise stated. This also implies faster performance since only the 
Zi- variables, with simplified updating equations, are needed. 



The complete algorithm will look as follows: 



1. Initialize z^ (and a.^ if present) randomly, suitably in the neighborhood of a truncated high- 
or low-T series solution. 

2. For each i: 
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• Update Zi (and a^) according to eqs. ( |39| ) with suitable step-sizes. 

• Correct all according to eq. (po|). 

3. Check if converged; if not, go to 2. 

4. Extract and Gij = z,; • Zj, and compute variational averages of interest. 

Typical step-sizes are k, 1/6 and ea ~ 1/2. The convergence check is done based on the rate 
of change and on the virial identity. The number of computations in each iteration step for this 
procedure is proportional to . The number of iterations required for convergence is a slowly 
growing function g{N). In total, thus, the execution time of the algorithm grows with N as N^g{N). 
In terms of CPU requirement convergence of a iV = 40 chain (a.; = 0) requires 3 seconds on a DEC 
Alpha workstation. 



4 High and Low T Results — Analytical Considerations 

Before embarking on a numerical evaluation of the variational approach with comparisons to MC 
results, it is interesting to see what can be gained from studying the high and low T limits, where 
analytical methods can be used. 

At the energy minimum, prevailing at T = 0, the polyelectrolyte will form a straight line. When the 
temperature is increased there will be a competition between the entropy and the repulsive Coulomb 
forces, and as T — > oo, the chain becomes Brownian, and the elongated or ordered structure is gone 
altogether. The temperature range of the transition from ordered to disordered structure is iV- 
dependent. As N increases at fixed T, the Coulomb force becomes relatively more important and 
the system effectively behaves as if T decreased: the polymer configuration becomes increasingly 
aligned. 

In the variational approach, as discussed above, the high temperature regime is characterized by a 
purely fluctuating solution, reflecting the Brownian nature of the chain at T — > oo. Such a solution 
survives as a local minimum also at lower T, where however also a rigid solution exists. Below a 
certain critical temperature Tc, the latter gives the global minimum, indicating a flrst order phase 
transition. This is probably an artefact of the variational approach - in the MC simulations the 
system shows no evidence of possessing a phase transition (see section 5). The rigid solution mirrors 
the ordered elongated structure of the polymer at low T. 

With these qualitative arguments in mind, we turn to a more detailed investigation of the behaviour 
of the polyelectrolyte in the high and low T limits, together with an evaluation of the corresponding 
variational results. The unscreened and screened cases will be treated separately. 



4.1 The Unscreened Coulomb case 
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4.1.1 High Temperature 



In the high T limit, the variational results can be expanded in l/T (see Appendix D). Thus, for the 
expectation value of the Gaussian energy Eq, the first two terms of the expansion yield, 

N-l _ 

(-Eg) =3(7V-l)r/2 + l/V2^ V (41) 

fc=i 

which agrees with the exact result to the order shown; the first discrepancy occurs in the 0{T~^) 
term, as is in fact true for any quadratic expectation value (riTj). By the virial identity, this also 
holds for the Coulomb energy (Ec) and for the total energy as well. 

For the individual bond lengths, the high T result can be written 

(r^) w ST + 4 V2/7rT (Vi + - //v) (42) 

where the last term is obtained from a continuum approximation, valid for large N. 



4.1.2 Low Temperature 

In the low T limit (see Appendix D), the exact result for the total internal energy, expanded in 
powers of T, is, 

(E) = Eo + {3N ~5)T/2 + 0{T^) (43) 

where Eq is the minimum energy at T = and 3 A'^ — 5 is the number of degrees of freedom, modified 
for the spherical symmetry (fi — 2). It turns out that with the above first order low T correction, 
the energy, and thus also Eq, Ec and (see below), are quite well approximated for the sizes 
and temperatures considered in this paper. 



At low T, the variational free energy is minimized by the rigid solution, for which the corresponding 
expansion is, 

{E)v = Eo + 3{N -l)T/2 + 0{T^) (44) 

where the first term is the same as in eq. (p3[), while the second term is qualitatively correct for 
large A^. 



Yet another low temperature expansion results from the purely fluctuating variational solution and 
it gives, 

{E)v = (6/7r)i/3£;o + 3(A^ - 2)T/2 + 0{T^) (45) 

which shows a 24% discrepancy in the first term. The same factor, (G/tt)^/'^, results for any quadratic 
expectation value and for any single term in (Ec) in this limit. For r.m.s. distances, like the 
monomer- monomer distance rmm. 



(46) 
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or the end-end distance Tee, 

this corresponds to an error of 11%. 

With the variational approach thus satisfactory in both temperature hmits one can hope to find it 
a reasonable approximation also at finite temperatures, as will indeed be borne out in section 5. 



4.1.3 Zero Temperature 



Having low T expansions under control in terms of the T = configurational properties, the latter 
remain to be calculated. They are given by the minimum energy configuration, which is aligned, 

Ti = bi-n (48) 

and unique (up to global translations and rotations). 

The bond lengths bi> satisfy the equation 

^» = E ^ (49) 

where b^ = Xljeo- length of the subchain a. This equation cannot be solved analytically 

(except for very small N), but a fair large- A'' approximation can be obtained (see Appendix E for 
details). This gives for a distinct monomer-monomer bond the result 

/ 2a/2 ^, 

Vi/T^O — ~ 

As a consequence, we find that at zero temperature the average bond-length should scale logarith- 
mically with N, 

rmm oc (logN)^/^ (51) 

For the end-end separation this implies 

reeOc7V(log7V)V3 (52) 

This result is interesting, since it predicts a scaling faster than N; the extra logarithmic factor comes 
from the stretching of the harmonic bonds. 

The results above might seem as rather academic and of little practical impact. However, the 
low temperature expressions turn out to be quite accurate when compared with MC results. In 
other words ordinary room temperature and aqueous solution corresponds to a surprisingly "low" 
temperature. This indicates that the low temperature expansion might be a good starting point for 
further work in polyelectrolyte solutions. 



i{N-i) 



log ( const 

V N 



(50) 
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4.2 The Screened Coulomb case 



At high T, an analysis similar to the one carried out for the pure Coulomb chain, can be done 
for the screened chain (see Appendix D). Also there, the variational approximations to quadratic 
expectation values turn out correct to next-to- leading order; the same holds for E and Eq, while 
Ec (which is one order down) is correct to leading order. 

Also at low T, the results remain essentially the same for the screened chain. Thus, the rigid solution 
gives correct energies to lowest order in T, while the purely fluctuating solution does not. 



5 Monte Carlo Simulation Techniques 



For the numerical evaluation of the variational approach, the results were compared to those from 
MC simulations, which were performed in the canonical ensemble with the traditional Metropolis 
algorithm | |3^ . For short chains this is a straightforward procedure, but for linear chains consisting of 
more than about 100 monomers convergence problems appear. Typically for a chain of 40 monomers 
four million moves/monomer were required in order for the statistical fluctuations in the end-end 
separation to be less than one per cent. The energy terms and local conformational properties 
like monomer-monomer separations converged much faster. The addition of salt, i.e. the use of a 
screened potential, improves the convergence characteristics, but on the other hand its evaluation 
is more time-consuming than the pure inverse square root. Careful coding, with table look-ups for 
the inverse square root routine and, in particular, for the screened Coulomb potential, turned out 
to be more rewarding, reducing the computation time by almost a factor of three. 



In order to treat longer chains with reasonable statistics we are forced to use more efficient algorithms 
like the pivot algorithm, first described in refs. jsj, with a high efficiency for linear chains on a 



lattice with short range interactions. Recently it has also been used successfully |36| for off-lattice 
simulations of a single polymer chain. It could be argued that the pivot algorithm should be even 
more efficient for chains with long range repulsive interactions like a charged polymer. The form 
of the pivot procedure used in this work can be described as a two step process, consisting of a 
random translation followed by a random rotation or vice versa. For a polymer chain with fixed 
bond lengths only random rotations will be used. The procedure is as follows: choose a monomer i 
and apply the same random translation to monomers j -I- 1 to TV. Then choose an axis at random and 
perform a random rotation of monomers i + l to N around this axis. Evaluate the interaction energy 
between monomers 1 to i, and i + 1 to A^. This is a quadratic process in contrast to the single move 
algorithm, which only requires the evaluation of N pair interactions/move. Finally a Metropolis 
energy criterion is used to test for rejection or acceptance of the new configuration. We find that with 
a maximal random displacement of the order of 5-10 rA and a maximal random rotation of tt, we 
reject approximately 50 % of the attempted moves. Typically, we generate 10^-10'* passes (one pass 
= one attempted move/monomer) resulting in a statistical uncertainty in the end-end separation of 
approximately one per cent. The uncertainty in the average monomer- monomer separation and in 
the Coulomb and Gaussian energies is much less. Local averages, however, like the ith bond length, 
may have larger uncertainties, something that is also discussed by Madras and Sokal |Q. The 
pivot algorithm seems to be superior to the traditional single monomer procedure described initially 
already for chains with A^ > 20. We also have found that restricting the procedure to translational 
moves still makes it superior to the traditional algorithm. The pivot algorithm makes it feasible to 
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simulate chains with more than one thousand interacting monomers. The major drawback with the 
pivot algorithm seems to be its limitations to linear chains or at least chains with simple topologies. 



Without excessive fine tuning of the translational and rotational displacement parameters, we find 
that the computational cost grows as N^. This power results from iV^ for each sweep of monomer 
moves, and an additional factor N from autocorrelations in quantities like end-end separations. 

Another efficient algorithm has recently been developed in ref. po|| . By identifying the slow modes 
in a Fourier analysis, one is able to use different random step lengths for different modes. This 
technique seems to be as efficient as the pivot algorithm and we have used it to check the accuracy 
of our simulations. For all cases investigated we obtain, within the statistical uncertainties, identical 
averages. The same is true for shorter chains where we also can use the original single monomer 
algorithm as a further test. 



6 Numerical Results 



The superiority of the purely fluctuating variational solution over the rigid one, as discussed in 
section 5, is illustrated in figs. and ||, where the variational results for ree and rmm are compared 
to MC data for iV = 20 and 80. 

The remainder of this section will be devoted to comparisons of the variational approach to MC sim- 
ulation results focusing on (1) energies, in particular the free energy and (2) various configurational 
measures. The = variational solution will be consistently used. 



6.1 Free and Internal Energies 

In table |]the variational results for internal Coulombic and Gaussian energies are compared with the 
MC results; the relative deviations are seen to increase both with increasing Cs and with increasing 
N. For fixed Cs the deviations seem to converge to constant values at large N. Hence it should be 
possible to extract "asymptotic" correction factors from comparisons at moderate N, do a variational 
calculation for a very large N, and predict what an MC calculation would give. 

A strong advantage of the variational approach is the direct access to the free energy, which is much 
more difficult to obtain in a MC simulation, requiring a cumbersome integration procedure. In 
order to evaluate the variational results we nevertheless attempt to estimate F(T) from MC data 
for iV = 20 using the following procedure. In dimcnsionlcss units one has 



Thus, we can define an excess free energy with respect to some reference temperature Tr as 




(53) 




(54) 
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Figure 1: (a) r^e and (b) as functions of T for an unscreened chain with N = 20. Filled 

circles represent MC data, and solid and dashed lines variational results, with a ^ and a = 0, 
respectively. 



which is then accessible in MC by a temperature integration of (E). In fig. |3|the excess free energy is 
shown as a function of T, for an TV = 20 chain with Tr corresponding to 1422 K. As can be seen the 
variational solutions for F{T) reproduce the extracted MC values very well in a wide temperature 
interval. Comparisons for larger iV are not feasible due to the indirect cumbersome MC extraction 
procedure. 

We remark that one possibly efficient alternative route to obtain free energies for different degrees 
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Figure 2: (a) r^e and (b) rmm as functions of T for an unscreened chain with = 80. Same 
notation as in fig. ||. 



of screening at fixed T in MC, would be to perform an integration in the Debye screening length and 
calculate the incremental excess free energy for a change An. The advantage of such a procedure 
would be that every point along the integration path corresponds to a physically realistic situation, 
while the completely screened chain would serve as a reference state. 
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N 


20 


40 


80 


160 


320 


512 


c,=0.0 Ec (V) 
(MC) 


6.20 
5.25 


7.58 
6.30 


8.80 
7.28 


9.94 
8.16 


11.0 
8.99 


11.7 
9.53 


Eg (V) 
(MC) 


6.65 
6.25 


7.40 
6.78 


8.08 
7.31 


8.66 
7.79 


9.20 
8.23 


9.54 
8.47 


c,=0.01 Ec (V) 
(MC) 


3.55 
2.70 


3.80 
2.83 


3.95 
2.89 


4.02 
2.91 


4.05 
2.93 


4.07 
2.93 


Eg (V) 
(MC) 


6.20 
5.70 


6.63 
6.00 


6.88 
6.15 


7.00 
6.22 


7.07 
6.26 


7.10 
6.27 


c,=0.1 Eg (V) 
(MC) 


1.90 
1.35 


2.03 
1.38 


2.16 
1.40 


2.19 
1.41 


2.15 
1.42 


2.15 
1.42 


(V) 
(MC) 


5.40 
5.00 


5.68 
5.15 


5.86 
5.26 


5.94 
5.29 


5.93 
5.32 


5.94 
5.32 


c,=1.0 £;c (V) 
(MC) 


0.65 
0.40 


0.70 
0.40 


0.74 
0.41 


0.75 
0.43 


0.76 
0.42 


0.76 
0.42 


Eg (V) 
(MC) 


4.35 
4.10 


4.53 
4.25 


4.61 
4.29 


4.67 
4.33 


4.69 
4.37 


4.70 
4.37 



Table 1: Average internal Coulombic and Gaussian energies per monomer in kJ/mol ■ monomer 
for unscreened and screened Coulomb potential. MC and V stands for Monte Carlo and variational 
calculations respectively. The salt concentration, c^, is given in molar (M). 

6.2 The End-End Separation 

The end-end separation is a critical measure of the accuracy, being a global quantity with contribu- 
tions from all bond-bond correlations. Unfortunately, it also turns out to be a complicated quantity 
from the convergence point of view in standard MC simulations, which means that it will have larger 
uncertainties than for example the average monomer- monomer separation - this is particularly true 
for the pure Coulomb chain. Table ^ contains a comparison between variational and simulation 
results for and r^e using the unscreened Coulomb potential. The result from the variational 
approach is impressive - the maximal deviation, 7.8% for N — 1024, from the MC results is well be- 
low the 11% bound discussed above. (Due to memory limitations on the local workstation N = 2048 
has not been pursued with the variational approach.) This result can be compared to the results 
from a (spherical) mean field approach [Q, which for the same system shows a deviation of 20% 
already for N = 100. In section 4, r^e/N and rmm were conjectured to vary linearly with (logiV)^/^ 
for large iV at zero temperature (cf. eqs. ( ^ , |2|)). In fig. ^, this hnear dependence is tested on 
both MC and variational data at room temperature (298K), with a surprisingly good result. This 
indicates that room temperatures can be considered low for a reasonably long Coulomb chain with 
the chosen parameters. It is also clearly seen how the variational approximation to r^e exceeds the 
MC values; asymptotically we expect an 11% discrepancy, as discussed above. 

The relative errors for Tf^^/N and r„im in table ^ may be used for a numerical estimate of the 
asymptotic error. Assuming that the relative error decays to the final value like {[ogN)~^/^ we 
estimate the asymptotic errors in r^e and r„im to be 10 and 13 %, respectively, in good agreement 
with the expected 11 %. 
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Figure 3: The excess free energy AF{T) (see eq. (p4[)) from the MC data (fiUed circles) and from 
the variational approach (solid line with a ^ and dotted line with a = 0) for (a) Cg = O.OM and 
(b) Cs = O.IM respectively (A^=10). 



The difference in slope for the variational and simulated Vee/N in fig. |jb is a consequence of the 
finite number of monomers. Numerically we find, at the present temperature, that N = 1000-2000 
is not enough in order to reach the asymptotic regime - the slopes in fig. ^ are larger than 1/3. 
By lowering T, however, we expect both slopes to approach this limiting value. 

The conjectured zero-temperature scaling results seem to contradict the results by deGennes et al. 
p8| and Baumgartner ||2^ that Vee should scale linearly with TV for an unscreened polyelectrolyte. 
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N 


20 


40 


80 


160 


320 


512 


1024 


2048 




13.04 


13.60 


14.11 


14.57 


14.99 


15.26 


15.63 




MC 


12.56 


13.01 


13.43 


13.81 


14.17 


14.38 


14.68 


14.99 


diff. 


3.8% 


4.5% 


5.1% 


5.5% 


5.8% 


6.1% 


6.5% 




Tee V 


122 


277 


632 


1425 


3152 


5340 


11478 




MC 


119 


269 


606 


1347 


2958 


4985 


10651 


22507 


diff. 


2.5% 


3.9% 


4.3% 


5.8% 


6.6% 


7.1% 


7.8% 





Table 2: and r^e in rA for unscreened Coufomb potential (c^ = O.OM) as computed with the 

variational (V) and Monte Carlo (MC) methods. The errors originating from the MC runs are 
estimated to be 0(0.2%). 



However, the latter result is valid only if the monomer-monomer bonds are rigid, but with elastic 
bonds as in the present model, eq. (|]), swelling is possible and ree/N increases slowly with N due 
to the long-range Coulomb repulsion. 

Returning to the simulations by Baumgartner ]20| , we note that his effective temperature is almost 
a factor of ten larger than in the present study. Such a high temperature means that the chain 
behaviour is essentially brownian in character making it numerically difficult to detect the electro- 
static expansion of the chain in a traditional MC simulation. The rigid monomer-monomer bonds 
used by Baumgartner also precludes the extra expansion predicted by eqs. (^ 52). 



Table I contains the same quantities as table ^ but for a screened Coulomb potential. The agreement 
detoriates, when a small amount of salt is added. The screening reduces the Coulomb repulsion, 
which in a sense is the hard part in our variational calculation, and one would naively expect an 
improved accuracy with a decreased interaction. The opposite result is found and the discrepancy 
in the end-end separation becomes as large as about 40% with 10 mM of salt and N = 160. At 
sufficiently high salt concentration the agreement improves again, as it should, with the Coulomb 
repulsion completely screened. With a strong screening the Coulomb potential will have a short 
range, and for sufficiently large N the chain will be Brownian. This is also reflected in tabled, where 
the agreement for r^e is worst for intermediate chain lengths and improve again when N increases. 
The large discrepancy seen e.g. for Cs=0.01 M and N—160 might seem surprising, considering the 
excellent agreement found for the pure Coulomb chain, but it reflects the difficulty to properly 
emulate a short-range potential with a harmonic effective energy |Tq] . 



6.3 Monomer-Monomer Separations 



The variational results for the average monomer- monomer separation is in excellent agreement 
with MC results for both the unscreened and screened cases - the largest error seen is of the order of 
5% (see tables ^and^). As for Vee the variational estimate is always larger than the MC value. This 
is also true for any single monomer-monomer separation (r^)^/^ as can be seen in fig. |^. The shape 
of the curves are correctly reproduced by the variational solution and the largest discrepancy is not 
unexpectedly found in the middle of the chain, which may be explained by a stronger accumulated 
electrostatic repulsion there. These results can be compared to the mean field solution of ref. p5[ , 
in which gradually more and more pair interactions were treated explicitly and withdrawn from the 
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1.4 



1.6 

[log(N)] 



1.8 

1 / 3 



1.05 



0.95 - 



0.75 

























• 







0.1 



0.25 0.4 
log (log N) 



0.55 



Figure 4: (a) as a function of (log A^)-'^/'^. Filled circles represent MC data and solid line 

the variational results. The dashed line is a linear fit to the MC data, (b) logree as a function 
of log(logiV). Filled and open circles represent MC data and variational results respectively. The 
lines are linear fits. 



mean field. It was found that only after the inclusion of next-nearest neighbours, leading to lengthy 
numerical calculations, the curve shape became qualitatively correct. 

We discussed above the zero-temperature scaling for an unscreened chain, r.,„„j cx 

(log AT) 1/3, When 

salt is added Tram increases much more slowly with chain length and for Cs=l M, when Tram 3> kT^ ■, 
it is essentially independent of chain length. One also notes that the individual bond lengths {rf)^^'^ 
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N 


20 


40 


80 


160 


320 


512 


Cs— U.Ul 




V 


12.60 


12.87 


13.02 


13.10 


13.14 


13.16 






A/rr 


12.09 


12.24 


12.31 


12.34 


12.36 


12.37 






Qin. 


4.2% 


5.1% 


5.8% 


6.1% 


6.3% 


6.4% 






V 


104 


201 


377 


680 


1 1 88 


1 71 






MC 




1 83 


31 7 


521 


895 


1110 






diff. 


4.4% 


9.8% 


19% 


31% 


44% 


54% 


c.=0.1 




V 


11.77 


11.90 


11.97 


12.01 


12.04 


12.04 






MC 


11.30 


11.35 


11.39 


11.38 


11.40 


11.40 






difF. 


4.2% 


4.8% 


5.1% 


5.5% 


5.6% 


5.6% 




^ee 


V 


78.2 


136 


231 


387 


640 


895 






MC 


72.9 


120 


192 


301 


459 


622 






difF. 


7.3% 


13% 


20% 


29% 


39% 


44% 


c.=1.0 


'^rnin 


V 


10.57 


10.69 


10.69 


10.69 


10.70 


10.70 






MC 


10.27 


10.29 


10.29 


10.30 


10.34 


10.32 






diff. 


2.9% 


3.9% 


3.9% 


3.8% 


3.5% 


3.7% 






V 


55.0 


86.9 


137 


217 


343 


468 






MC 


52.1 


79.5 


122 


182 


283 


364 






diff. 


5.6% 


9.3% 


12% 


19% 


21% 


29% 



Table 3: rmm and r^e in rA for the screened Coulomb potential (cs=0.01M,Cs=0.1M and Cs=1.0M 
respectively) as computed with the variational (V) and Monte Carlo (MC) methods. The errors 
originating from the MC runs are estimated to be 0(0.1%) 



do not vary in the central part of the chain when salt is present (see fig. |5|). 

Furthermore, for an unscreened chain, following eq. (|5C|), we expect at T = the individual bond- 
lengths, raised to the power three, to vary linearly with u = log[s(l — s)], where s = i/N. In 
fig. I^a it is shown that this is approximately true, more so for the variational solution than for 
the MC results. The lower curve in fig. obtained with a screened Coulomb potential, shows a 
qualitatively different behaviour. Fig. |6|3 contains a similar graph for different TV. 

The T — scaling relation for individual bonds, eq. (|50|), has the peculiar consequence that the 
length of a bond at the end of the chain becomes independent of N. This can be seen by rewriting 
the scaling relation as 

(r?)i/2 « {logN + log[s(l - s)]}i/3 « [logz]i/3 (55) 

where the last expression holds for small i. Eq. ( p5| ) also holds for the MC results, where we find 
the first few bond lengths to be independent of N. 
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10 20 30 40 

Bond (i) 



Figure 5: Bond lengths (rf)-'^/^ along an = 40 chain. Salt concentrations Cg from top to bottom 
are M, O.OIM, O.IM and 1 M respectively. Solid line represents variational results and dashed line 
Monte Carlo data. 



6.4 Angular Correlations 

In order to further test the variational solutions, we also have calculated angular correlations between 
bonds, 



which roughly gives the average of the cosine between bonds. In fig. |^a, variational and MC data 
is shown for the neighbor correlation Ci^i-^-i. It is seen that the variational Ansatz consistently 
overestimates the angle, more so in the presence of salt than without. Comparing figs. |^ and|^a, we 
find that the above discussed discrepancy between the MC and variational results for the end-end 
separation seems to be due to differences in angular correlations as well as in bond lengths. For the 
unscreened case the two sources seem to be of comparable magnitude, while for the screened case 
the angular correlations seem to be larger and the main cause of discrepancy. 

Fig. I^b shows a more global angular correlation Ci^i between the first and all successive bonds. The 
variational results for this quantity are in excellent agreement with MC data, both in the screened 
and in the unscreened chain. 




(56) 
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- 4 - 3 - 2 - 1 

log[s(l-s)] 




log[s(l-s)] 

Figure 6: The cubed bond length (rf)'^/^ as a function of log(s(l — s)), with s = i/N the relative 
monomer position, (a) Upper line shows variational results and symbols MC data for N = 40, no 
screening. The zero temperature scaling corresponds to a straight line. The lower line shows MC 
data for a screened chain (c^ = 0.01 M) for comparison, (b) Variational results for unscreened 
chains of varying size N. 

7 Summary and Outlook 

A deterministic variational scheme for discrete representations of polymer chains has been presented, 
where the true bond and Coulombic potentials are approximated with a trial isotropic harmonic 
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Bond (i) 



0.2 




Bond (i) 

Figure 7: Angular correlations (a) between neighbouring bonds, and (b) between the first and 
successive bonds. Solid lines represent variational and dashed lines MC results. The upper pair of 
curves are for an unscreened chain, while the lower pair is for Cg = 0.01 M. {N = 80.) 



energy. The variational parameters obey matrix equations, for which a very effective iterative 
solution scheme has been developed - the computational demand is N^. 

The high and low T properties of the variational approach has been analyzed with encouraging 
results. Also, the approach is shown to obey the relevant virial identities. 

In contrast to MC simulations, the free energy is directly accessible with the variational method. 
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When confronting the results from the method with those from MC simulations, very good agreement 
is found for configurational quantities in the case of an unscreened Coulomb interaction (the error 
is within 11 %). 

In the screened case the method does not reproduce the MC results equally well although the 
qualitative picture of conformational properties is there. We attribute this problem to the difficulty 
for a Gaussian to emulate short range interactions. 

Recently, MC simulations were pursued for titrating Coulomb chains ||39|| . For such systems the 
Coulomb potential of eq. (Ff) is modified to 



where the binary variables Si are either 1 or depending whether monomer i is charged or not. 
Thus minimizing E now also includes a combinatorial problem - deciding where the charges should 
be located. Variational techniques related to the ones used in this paper have been successfully 
used in pure combinatorial optimization problems ||4C| ], where again tedious stochastic procedures 
are replaced by a set of deterministic equations. Along similar lines, the approach of this paper can 
be modified to allow for a variational treatment also of the titrating problem. 

The variational approach is also directly applicable to more general topologies - bifurcations pose 
no problems. Proteins could also be treated in this way provided the traditional Lcnnard- Jones 
potentials are replaced by forms that are less singular at the origin. 
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Appendix A. The Variational Approach — Generalities 



In this appendix wc discuss the generic variational approach that is used in this paper for the 
particular problem of a Coulomb chain. We consider a generic system, with dynamical variables x 
in some multi-dimensional state-space, and assume that a real energy function E{x) is given. 

For an arbitrary probability distribution P{x) in an arbitrary state space, the free energy with 
respect to an energy E{x) is generally defined as 

F={E)- TS (Al) 

where S is the entropy, 

S=-{\ogP) (A2) 
and expectation values are defined with respect to P. Writing P{x) as 



P{x) = ^eM-Ev{x)/T) (A3) 



with 



Z = j dx exp{-Ev{x)/T) (A4) 
the free energy can be written as 

F = -TlogZ + {E - Ev) (A5) 



Note that 

\ T 



exp(-F/r) = Zexp/ ^^^ ^ ) (A6) 



= exp{-F/T)\E^=E 

where the inequality is due to the convexity of the exponential function. Thus, F is bounded 
from below by its value for Ey = E, corresponding to the proper Boltzmann distribution, P{x) oc 
exp{-E{x)/T). 

The variation of F due to a variation SEy is given by 

SF = -T{{SEv{E- Ev)) - {SEv){E- Ev)) (A7) 
= -T{5Ev{E-Ev))c 
where {ah)c stands for the connected expectation value (cumulant), {ah) — {a){b). 

The idea of the variational approach is to choose a suitable simple Ansatz for the variational energy 
Ev, with a set of adjustable parameters a^, i = 1, . . . , Np, the values of which are to be chosen so as 
to minimize the variational free energy. Demanding the vanishing of the variation of the free energy 
due to variations in the parameters aj then leads to the general equations for an extremum: 

dEv ' ^ 



{E-Ev)) =0, i = l,...,Nj, (A8) 



c 
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where 0^ denotes an expectation value based on the variational Boltzmann distribution. This 
determines the optimal values of the parameters. Exact expectation values are then approximated 
by the corresponding variational ones. 
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Appendix B. The Variational Free Energy for the Chain 



In this appendix we derive the expressions for the variational free energy (eqs. (|2^, ^)), for the 
unscreened as well as the screened Coulomb chain, with or without translational parameters in the 
variational Ansatz. 



Unscreened Coulomb Chain 

For the specific case of the Coulomb chain of length N, the energy amounts to 

i a- 

where cr is a contiguous subchain. 
Gaussian Parameters Only 

We consider first a pure Gaussian variational Boltzmann distribution, corresponding to 

hi 

The parameter matrix G is forced to be symmetric and positive-definite by expressing it as 

The general expression for the variational free energy is 

F ^ -TlogZv ~ {Ev)v + {E)v (B4) 

where the expectation values are with respect to the normalized variational Boltzmann distribution 
exp(— i?y/r)/Zy. The first two terms are trivial to compute. The first is 

ST 

-TlogZv = — logdetG^i = -3Tlogdetz (B5) 
apart from a trivial constant that can be neglected, as can the second term, 

- {Ev}v ^ ^1{N - 1)T (B6) 

The last term, 
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consists in a sum of terms, each amounting to the variational expectation value of a simple function 
of a Gaussian vector variable To-, of which is a special case. Its probability distribution is given 

by 

P(r.)(xexp(^-^^ (B8) 

with 

= ^ z, (B9) 

Thus, we have 

{r-)v = 3z2 (BIO) 

and 

Summing up, the variational free energy takes the form 

F = -STlogdet ^ + ^ z2 + ^ U^{z^) (B12) 

i <y 

to be minimized with respect to the variational parameters z^. 
Gaussian and Translational Parameters 

For the more general variational Ansatz with additional translational parameters aj , 

the main difference is a translation of the individual probability distributions of eq. ( [B§| ), which 
now read 

J^(r.)(xexp(^- ^^-~J-^' ^ (B14) 

with 

a, = ^ a, (B15) 

This gives 

(r2)v = 3z2 + a2 (B16) 

and 

1\ ^^erff-^) ^[/f(z.,a.) (B17) 



The variational free energy becomes 



F = -STlogdet z + i ^{3z^ + af) + Y, U^{za,a^) (B18) 
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Screened Coulomb Chain 



Next we consider the Debye-screened version of the Coulumb chain, with the energy 

^ l^^,^^exp(_«^ (B19) 



Gaussian Parameters Only 



For the variational Ansatz with only Gaussian parameters, eq. ( ]B2| ), we need eq. (BIO) and the 
expectation value 



exp(-Krcr) 



'2 1 

= \l K CXp 



2 ^erfc^-|)^C/f(z.) 



(B20) 



The variational free energy then reads 



F = -3T log det z + - ^ z,^ + ^ C/f^ (z,) 



(B21) 



Gaussian and Translational Parameters 



Finally, if for the screened Coulomb chain, eq. (B19), also translational parameters are used in Ey, 
eq. ( |B13 ), we will need the following result in addition to eq. (BIE), 



/ exp(-KrCT) 
\ 



exp(At^zg/2) 

2acr 
U2{Za,acr) 



exp(— Kacr)erfc 



( 



— ) — exp(Kao-)crfc 



\/2zg 



The variational free energy will read 

F = -STlogdet Z + ]- ^(3z2 + a.^) + ^ U^{z,,a,) 



(B22) 



(B23) 
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Appendix C. The Virial Identity 



In this appendix we derive the virial identities and show that these are respected by the variational 
approach. 



Exact Virial Identity 



For any system described by a Boltzmann distribution 



P(x) = lexp(-i;(x)/T) (CI) 



with X e TZ^ , and E rising as a power for large |x|, we will have 



y" V • (f (x) exp(-£;/T))dx = 



(C2) 



Z 

for e.g. any polynomial f , due to the integrand being an exact divergence. This is equivalent to 

T(V • f) = (f • V£) (C3) 
Thus, by varying f , we can obtain an infinite set of identities for the system. 

The virial identity results from the particular choice f = x; in its general form it reads 

(x • VE) = TD (C4) 
where D is the dimension of x-space (x • V is the scaling operator) . 

This is particularly useful if the energy E is given by a sum of terms Ea homogeneous in x, 

X • = XaEa (C5) 

in which case the virial identity takes the simple form 

Y,K{Ea)=TD (C6) 



This applies e.g. to the case of the unscreened polyelectrolyte. There the scaling operator is given, 
in relative coordinates, by X^^rj • Vr^, and we have Ag = 2 and Ac = — 1; the virial identity thus 
reads 

2{Eg) - {Ec) = 3(iV - 1)T (C7) 



Vctriational Virial Identity 

The virial identity is preserved by the variational approach under certain conditions, to be specified 
below. For the generic system above, minimizing the free energy 

F = Fv + {E-Ev)v (C8) 
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w.r.t. the parameters a; of a variational energy Ey, leads to (cf. eq. (A8) in Appendix A) 



(C9) 



Now, choosing f — x(i? — Ey ) in cq. (03) with the variational Boltzmann distribution, we have 

DT{E - Ev)v + r(x • V(E - Ev))v = {{E - -By)x • \'Ev)v (CIO) 
On the other hand, since the virial theorem holds for Ey, 

DT = (^■VEv)v (Cll) 



Substituting this into eq. (CIO), we obtain 

r(x • ViE - Ev))v = {{E - Ev)^ ■ VEv)l 



(C12) 



If the set of parameters of Ey is such (and this is the crucial condition) , that the scaling operation 
on Ey can be written in terms of derivatives with respect to the parameters, i.e. if 



y^-VEv = ^Gj(c 



^dEv 
dai 



(C13) 



then the righthand side of eq. (C12) vanishes at the minimum, due to eq. (C9), and we are left 
with 

{^■VE)v = {yi-VEv)v = DT (C14) 

which is what we desired. 

Note that the derivation only relies on a local extremum of the free energy. Thus, for the polymer. 



the virial identity, eq. (C7), is respected by both the rigid (a ^ 0) and the purely fluctuating (a — 0) 
solutions. 
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Appendix D. High and Low T Expansions 



High T expansions 
Exact results 

At high T, the chain size wiU be large, and accordingly, the Gaussian term will dominate over the 
interaction V in the energy expression, 

E = EG + V=lj2^^+J2vir<r) (Dl) 

i a 

with a summed over contiguous subchains. 

It is then natural to attempt an expansion in the perturbation V. For an arbitrary expectation 
value, we have the perturbative expansion 

(/) = if)' - ^{fV)% + ^.{fVVfc (D2) 

where ( )^ refers to connected expectation values (cumulants) in the unperturbed Boltzmann distri- 
bution. Due to the singular behaviour of v{r) for small r in the (screened or unscreened) Coulomb 
case, only the first few terms will be finite. This indicates that expectation values cannot be ex- 
panded in a pure power series in T, and that logarithmic corrections will occur after the first finite 
terms. 

We are interested in quadratic expectation values of the type (r^ • r^). These can be combined to 
give e.g. the rms cnd-to-ond distance Vee, the gyration radius, and the Gaussian energy {Eg) (and 
thereby, in the pure Coulomb case, the interaction energy (Ec) by the virial identity). 

For the pure Coulomb chain, v{r) = 1/r, we get the perturbative expansion 

where a denotes a contiguous sub-chain containing the iih and the jth bond, and io- its total 
number of bonds. This leads to 

{Eg) = 3{N - l)T/2 + ^lil^ 2.-1/2 ^ ^(T-^) (D4) 

(T 

and 

rl = 3(iV - 1)T + ^ ^1/2 + o(r-2) 

(T 

Similar results are obtained for the case of a screened Coulomb interaction, v{r) = e~^^/r, where 
the expansion of a quadratic expectation value gives 

(r, . r,) = 3T% + -^^A ^ ^-5/2 + o{T-') (D6) 
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Variational results 



The corresponding variational results can also be expanded at high T (where a^ = 0), by expanding 
the variational solution z^^ around the unperturbed value, which can be chosen as \/Tdifj_. The 
variational approximation to a quadratic expectation value, (r^ • rj)y = 3zi • Zj, can then easily be 
expanded. The results thus obtained reproduce the exact results, eqs. 
order shown. 

We conclude that, independently of screening, the high T variational results are correct to next-to- 
leading order for E and Eq, and to leading order for Ec- 



Low T expansions 

Here we will treat only the pure Coulomb case in detail; most of the discussion applies also to the 
screened case. 



D3.D6), correctly to the 



Exact results 



At low T, expectation values can be expanded around the configuration that minimizes the energy. 
This expansion is slightly complicated by the the rotational degeneracy of the minimum. The 
results can be expressed in terms of the classical configuration, which is given by a straight line 
configuration. 

Let bi be the the bond-lengths at the energy minimum. These have to be computed numerically, by 
solving the equation 

where b„ is the length of a subchain containing the ith bond. Note that the above equation can be 
written as a matrix equation: 

Thus, h is an eigenvector of the matrix B with a unit eigenvalue. Similarly, we can define a whole 
series of tensors: 



Ec 



Ef 

E^^^' (DIO) 



E i (DID 
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They are all symmetric, and contracting either with bi gives the tensor of rank one less. 



In addition, we need two more matrices, related to B, 

U = (l + 2B)-i (D13) 
V = P{1-B)-^P (D14) 

where P denotes the projection matrix onto the subspace orthogonal to b, which is deleted by 1 — i3. 
In terms of these tensors, we have the quadratic expectation- values at low T: 

Iim-Vi„,)]+OiT^) (D15) 



(r, -rj) = b,bj +T i U,j + 2V^j + + 3 ^ Ckl■m{b^UJk + b,U^k){Uh 



where the first two terms of the T coefficient are the naive contributions from the longitudinal and 
transverse fluctuations. The rest are corrections due to the rotational degeneracy of the T = 
configuration, which is also responsible for the transverse zero- modes (of 1 — B). 

From this we obtain e.g the average Gaussian energy, 

{Eg) = l/iEo + T{3N/2 - 11/6) + 0{T^)q (D16) 
from which, using the virial identity, we obtain 

(Ec) ^2/3Eo-2T/5 + 0{T^) (D17) 

and 

(E) ^ Eo + T{3N ~5)/2 + 0{T^) (D18) 

where Eq — 3/2^^6^ is the exact energy at T = 0. In the last equation, the T-coefficient is, as it 
should, half the number of degrees of freedom, not counting the two rotational zero-modes (and the 
three translational ones already removed). 



In the screened case, similar results can be obtained. In particular, eq. (D18) remains valid, though 
with a different Eq. 



Variational results 



Similarly, the variational results can be expanded at low T. We have to distinguish between the 
two different solutions. 

For the purely fluctuating solution with a^ = 0, the variational free energy is, at T = 0, 

This is obviously just the energy of an (N — l)-dimensional version of the chain, with modified 
coefficients. It is minimized by the aligned configuration 




(D20) 
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where bi are given by eq. (D7), and n is some unit vector in iV — 1 dimensions. For small but finite 
T, we have to add the entropy term, --3Tlogdet{z}, to Fq. This forces the configuration out of 
alignment. The resulting configuration can be obtained as a low-T expansion around the T — 
solution. The first correction to z will be of order Vt, and since the T = are aligned, the 
matrix inverse diverges - it will go like 1 / \/T. 

For the total energy, the leading correction can be obtained as follows. The equation to solve is 

STw, = V,/'o(z) (D21) 

For the T = solution, ViFo(zo) — 0. Thus, to lowest order, 

3Tw, = V, ^ VjFo(zo) ■ dz, (D22) 
j 

The leading energy correction will be 

1 3T 
dFo = -Y, V,Vji^o(zo)c?Zjdz» — ^ w, • dz, (D23) 

ij i 

Now, dzi consists of an aligned part and a transverse part, both cx Vt, while for the aligned 
part is cx 1, while the transverse part is cx I/Vt. Thus, the leading contribution to dFo comes from 
the transverse part. Taking the trace of the tranverse part of the identity Wi • zj — Sij, and noting 
that the transverse part of z sits entirely in dz, we have to leading order 

d^o = ^^^^'^^ =d{E)v (D24) 



Using the virial identity, which holds also for the variational expectation values, we get to first order 
in T: 



0{T^) (D25) 



TT Z 

{Ea)v = l{-y/'Eo + ^^^^^ + 0{T^) (D26) 
on Z 

{Ec)v = l{-y/'Eo-T + 0{T^) (D27) 

3 TT 

where Eq is the exact T = energy. Note that already the zero-order results are off by a factor 
j-|.-ji/3 1.24, but that the correction to E is correct in the high N limit. 

Simila r res ults are obtained for the screened Coulomb chain. Most of the general analysis leading 
to eq.(D24) still holds, and we have e.g. 

{E)v^E', + ^^^^^^ + 0{T') (D28) 

with E'q^ Ea. 

For the symmetry-broken a; ^ solution, the T = configuration is instead given by 

a, = 6,n , z, = (D29) 
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with hi as in eq. (D7), and n an arbitrary unit vector in TZ^. The small T corrections are obtained 
from expressing F as 

F = -3riogdet{z} + (^(a, + ^ z,^3^))j (D30) 
where are uncorrelated standard Gaussian noise variables. Expanding in z, we obtain 

F = -3riogdet{z} + E{a) + ^i^' ' '^J^(^) + ■ • • (^31) 

Because the Coulomb term satisfies Laplace' equation, this is just (provided a 7^ 0) 

3 

F = -3T log det{z} + £'(a) + - ^ z,z, (D32) 

i 

and the variational free energy separates in z and a for small z. Thus, minimum is obtained for 

ai = nbi (D33) 



and z satisfies 

which means 

The energies become 



' STw, + 3zj = (D34) 
Zi ■ Zj = TSij (D35) 



{E)v = Eo + ^^^^ (D36) 
(EaW - 1^0^'-^^ (D37) 
{Ec)v = "^Eo (D38) 
which is correct to lowest order, and for large N qualitatively correct to first order (except for Ec). 

Again, the screened chain lead to similar results; in particular, the T = energies will be the correct 
ones. 
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Appendix E. Zero Temperature Scaling Properties 



The T = configuration of a pure Coulombic chain cannot be obtained analytically, but must be 
computed numerically. However, an approximate calculation can be done. The equation for the 
bond lengths bi in the elongated ground state configuration is given by eq. (D7): 



(El) 



By assuming that the bond length is locally approximately constant, this can be approximated by 



' k=i i=i 



This can be rewritten as 



N-l 



1=1 



This in turn can be approximated by an integral, leading to 



log const - 



N 



Defining s — i/N, this amounts to 



bi « (log(const Ns{l - s))) 



1/3 



(E2) 



(E3) 



(E4) 



(E5) 



Eqs. (gH) give a quite accurate picture of the variation of the bond lengths at T = 0. They also 
imply, that the typical ground state bond length should swell roughly as (logiV)^/'^ for large N. 
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